Numerical simulation of transmission coefficient using c-number Langevin equation 
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. We numericaUy implement the reactive flux formalism on the basis of a recently proposed c- 

' number Langevin equation [Barik et al, J. Chem. Phys. 119, 680 (2003); Banerjee et al, Phys. 

Rev. E 65, 021109 (2002)] to calculate transmission coefficient. The Kramers' turnover, the 
enhancement of the rate at low temperatures and other related features of temporal behaviour of 
the transmission coefficient over a range of temperature down to absolute zero, noise correlation 
{DX)| and friction are examined for a double well potential and compared with other known results. This 

^ . simple method is based on canonical quantization and Wigner quasiclassical phase space function 

' and takes care of quantum effects due to the system order by order. 

^ '. 

I ^- I. INTRODUCTION 

^ : 

The classic work of Kramers^ on the diffusion model of classical reaction laid the foundation of the modern dynamical 
theory of activated processes. Since then the field of chemical dynamics has grown in various directions to extend the 
Kramers' result to non-Markovian friction?^, to generalize to higher dimension^*^ and to semiclassical and quantum 
(— I , rate theoriesSiSi^iiSiiiiiSiiii^iiSiiSiiLiLiS to apply extensively to biological processesSS and to several related issues. 
^ ' An important endeavor in this direction is the formulation of reactive flux theorji24i2iSii2& which has been worked 
^ [ out in detail over the last two decades by several groups. The method is essentially based on the realization of rate 
O ' constant of a chemical reaction as a Green-Kubo correlation function calculated at the transition state which acts 
as a dividing surface between the reactant and the product states. In the spirit of classical theory of Kohen and 
Tannor--^ we have recently formulated a quantum phase space function approach to reactive flux theory to derive a 
' transmission coefficient in the intermediate to strong damping regime^^. The object of the present paper is threefold: 
, Ph , (i) to extend the treatment to numerical simulation for calculation of time-dependent transmission coefficient (ii) to 
analyze the Kramers' turnover, the enhancement of the rate at low temperature and some related features of 
temporal behaviour of the time-dependent transmission coefficient over a wide range of friction, noise correlation and 
^ ' temperature down to vacuum limit and (iii) to confirm the validity of our analytical resultiS in the spatial diffusion 
limited regime. 

The present scheme of simulation is based on our recently proposed c-number Langevin equationi^*iSiiLiSiiS coupled 
to a set of quantum dispersion equations which take care of quantum corrections order by order. In what follows 
, we make use of these equations to follow the evolution of an ensemble of c-number trajectories starting at the top 
■ of the barrier and sampling only the activated events and recrossing according to classical numerical reactive flux 

formulation of Straub and BernaS^iSSiSS. 
' The outlay of the paper is as follows: we first give an outline of c-number Langevin equation in Sec. II followed by 
O [ a discussion on the method of numerical simulation for calculation of transmission coefficient using c-number reactive 
■ ^ . flux formalism in Sec. III. The numerical results are compared with those of others and discussed in Sec. IV. The paper 
J>-j' is concluded in Sec.V. 

Ph. 
> 

X 

' To start with we consider a particle coupled to a medium consisting of a set of harmonic oscillators with frequency 
LOi. This is described by the following Hamiltonian. 
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II. A BRIEF OUTLINE OF C-NUMBER LANGEVIN EQUATION AND THE COUPLED QUANTUM 

DISPERSION EQUATIONS. 



(1) 
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Here the masses of the particle and the reservoir oscillators have been taken to be unity, q and p are the coordinate 
and momentum operators of the particle, respectively and the set {qj,pj} is the set of co-ordinate and momentum 
operators for the reservoir oscillators linearly coupled to the system through coupling constant kj. V{q) denotes the 
external potential field which, in general, is nonlinear. The coordinate and momentum operators follow the usual 
commutation relations [q,p\ — ih and [(jijPj] = ihSij. 

Eliminating the reservoir degrees of freedom in the usual way we obtain the operator Langevin equation for the 
particle 



q{t)+ / dt'j{t-t')qit') + V'{q) = F{t) 
Jo 



(2) 



where the noise operator F(t) and the memory kernel j{t) are given by 



F{t) = ^ [^{9i(0) — q{0)}Kj cosLUjt + pj{0)Ky^ sinujjt 



(3) 



and 



or in the continuum limit 



7(i) = Kj cosojjt 
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lit) 



k{uj)p{uj) cosujt duj 



(4) 



piuj) represents the density of the reservoir modes. 

Eq.Q is an exact operator Langevin equation for which the noise properties of F{t) can be defined using a suitable 
canonical initial distribution of bath co-ordinates and momenta. In what follows we proceed from this equation to 
derive a Langevin equation in c-numbers. The first step towards this goal is to carry out the quantum mechanical 
average of Eq. 10) . 



m + / dt'^it - t'm')) + {v'm = (f^) (s) 



where the average (....) is taken over the initial product separable quantum states of the particle and the bath oscil- 
lators at t = 0, |0){|ai)|a2) |aAr)}. Here \4>) denotes any arbitrary initial state of the particle and \ai) corresponds 

to the initial coherent state of the i-th bath oscillator. \ai) is given by \ai) = exp(— |aip/2) ^^^p(a"'/\/niT)|ni), 
ai being expressed in terms of the mean values of the co-ordinate and momentum of the i-th oscillator, (9i(0)) — 
{Vh/2uji){ai + a*) and (pi(0)) = iy/huji/2 {a* — at), respectively. Now {F{t)) is a classical-like noise term which, 
in general, is a nonzero number because of the quantum mechanical averaging over the co-ordinate and momentum 
operators of the bath oscillators with respect to initial coherent states and arbitrary initial state of the particle and 
is given by 

{F{t)) = [ii^M) - im)}^j cosu;,t + {p,{0))k'/' sincujt] (6) 
j 

Now the operator Langevin equation can be written as 



(q)+ / dt'^{t-t'mt')) + {V'iq))=fit) (7) 







where {F{t)) = fit), denotes the quantum mechanical mean value. 

We now turn to the ensemble averaging. To realize f(t) as an effective c- number noise we now assume that the 
momentum {pj{0)) and the co-ordinate ((7j(0)) — (g(0)) of the bath oscillators are distributed according to a canonical 
thermal Wigner distribution for shifted harmonic oscillator^, as. 



3 



SO that for any quantum mean value Oj{{pj{0)), {((?j(0)) — (q(0))}), the statistical average (....)s is 



(o,), - J o, V, d{p,{Q)) d{{q,m - mm (9) 

Here fij indicates the average thermal photon number of the j-th oscillator at temperature T and is given by 
Bose-Einstein distribution fij = 1 / [exp[huj j / kT) — 1] and M is the normalization constant. 

The distribution Eq.® and the definition of the statistical average over quantum mechanical mean values Eq.Q 
imply that f{t) must satisfy 



{f{t))s = (10) 



and 



{f{t)f{t'))s ^\Y.^1 ^= COS^j(t-<') 

j 



or in the continuum limit 



{f{t)f{t'))s = 2 y dio p{lo) hcj f coth ^ j cosu;{t - t') 

^ c{t-t') (11) 

That is, c-number noise f{t) is such that it is zero-centered and satisfies the standard fluctuation-dissipation relation 
as known in the literature. For other details we refer to [16-20]. 

We now add the force term V'{{q)) on both sides of the Eq.© and rearrange it to obtain 

m + V'{q) + /* dt' 7(i - t') qit') = f{t) + Q{t) (12) 



where we put {q{t)) — q{t) and {q{t)) = p{t) ; q{t) and p{t) being quantum mechanical mean values and also 



Q{t)^V'{q)-{V'{q)) (13) 
represents the quantum correction to classical potential. 

Ea. H12|) which is based on the ansatz Eq.©, is governed by a c-number noise f{t) due to the heat bath, characterized 
by Eg. ljlOII . Ea. (|ll|l and a quantum correction term Q{t) characteristic of the nonlinearity of the potential. The 
canonical thermal Wigner distribution - the ansatz Eq.® is always positive definite. It goes over to a pure state 
distribution for the ground state of shifted harmonic oscillator in the vacuum limit, i.e. , at T = which is again 
a well-behaved^l distribution. The quantum nature of the dynamics thus arises from two sources. The first one is 
due to the quantum heat bath. The second one is because of nonlinearity of the system potential as embedded in 
Q{t) of Ea. (|15|) . To make the description complete for practical calculation one has to have a recipe for calculation 
of Q(i^)iSiiSiiLi24i2i2Si2i. For the present purpose we summarize it as follows: 

Referring to the quantum mechanics of the system in the Heisenberg picture one may write. 



q{t) = {q{t))+m) 
Pit) = {pit))+6pit) 



(14) 
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Sq(t) and Sp(t) are the operators signifying quantum corrections around the corresponding quantum mechanical 
mean values q and p. By construction {Sq) — {6p) — and [5q,Sp] — ifi. Using Eq.lQ^ in {V'{q)) and a Taylor series 
expansion around {q) it is possible to express Q{t) as 



QW = -E;^^^"^'^(9)('5'?"W> (15) 

n>2 

Here is the n-th derivative of the potential V{q). To second order Q{t) is given by Q{t) — ~^V"'{q)(5q'^{t)) 

where q{t) and {dq'^{t)) can be obtained as explicit functions of time by solving following set of approximate coupled 
equations Eq.((Tn|) to Eq.ljTHl together with Ea. lfH)) 

j^{6f) = (SqSp + SpSq) (16) 

j^{6qSp + 6p6q) - 2{6f) - 2V"{q){6f) (17) 

j^{6p^) = -V"{q){5q5p + 5p5q) (18) 

While the above set of equations provide analytic solutions containing lowest order quantum corrections, the 
successive improvement of Q{t) can be achieved by incorporating higher order contribution due to the potential V{q) 
and the dissipation effects on the quantum correction terms. In Appendix A we have derived the equations for quantum 
corrections upto fourth order— as employed in the present numerical scheme. Under very special circumstances, it has 
been possible to include quantum effects to all ordersiSi^. The appearance of Q{t) as a quantum correction term due 
to nonlinearity of the system potential and dependence of friction on frequency make the c-number Langevin equation 
fEa. (|12|l ') distinct from the earlier equations of Schmid and Eckern et aJ^24£. The approach has been recently utilized 
by us to derive the quantum analoguesi^*iSiiLi& of classical Kramers, Smoluchowski and diffusion equations with 
probability distribution functions of c-number variables. An important success of the scheme is that these equations 
of motion for probability distribution functions do not contain derivatives of the distribution functions higher than 
second for nonlinear potentials ensuring positive definiteness of the distribution functions. This is in contrast to usual 
situations'^* where one encounters higher derivatives of Wigner, Glauber-Sudarshan distribution functions for nonlinear 
potential and the positive definiteness is never guaranteed. More recently the classical Kohen-Tannor formalism of 
reactive flux has been extended to quantum domainiS. In what follows we present a numerical simulation of reactive 
flux using c-number Langevin dynamics. 



III. NUMERICAL SIMULATION 



The numerical solution of Ea. H12() . along with Eas. ljAl| to Ea. ljA3p is performed according to the following main 
steps. 

We first briefly outline the method of generation of c-number noise. Ea. (|ll|l is the fluctuation-dissipation relation 
and is the key element for generation of c-number noise. {f{t)f{t'))s is correlation function which is classical in form 
but quantum mechanical in content. We now show that c-number noise f{t) is generated as a superposition of several 
Ornstein-Ulhenbeck noise processes. It may be noted that in the continuum limit c{t — t') is given by 



c{t ^ t') = — I dio k{uj) p{uj) hu) I coth — - — j cosa;(i — t') (19) 
2 Jo V 2kT J 

In determining the evolution of stochastic dynamics governed by Eq. (|12|1 it is essential to know a priori the form of 
p{uj)k{lo). We assume a Lorentzian distribution of modes so that 



«:(^)p(c.) = - ( \ (20) 

where F and Tc are the dissipation in the Markovian limit and correlation time, respectively. Ea. H20() when used in 
Eq. in the continuum limit yields an exponential memory kernel 7(i) = (r/Tc)e-*/^'=. For a given set of parameters 
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r and Tc along with temperature T, we first numerically evaluate the integral (19) as a function of time. In the next 
step we numerically fit the correlation function with a superposition of several exponentials, 

c(t - = ^ — exp f ±^ V i = 1,2,3... (21) 

. Ti \ Ti J 

The set Di and the constant parameters are thus known. In Fig.l we compare the correlation c{t) determined 
from the relation (19) with the superposition (21) for three different temperatures at kT — 1.0, 0.5 and 0.1 for F = 1.0 
and Tc = 3.0. The numerical agreement between the two sets of curves based on Ea. (|19|l and Ea. H21|l suggests that 
one may generate a set of exponentially correlated color noise variables r]i according to 

r^^^-- + - m (22) 

Ti n 

where 

Mt))=Q and (6(0)ej(T)) =2A5y '5(t) (z = 1, 2, 3....) (23) 

in which £,i{t) is a Gaussian white noise obeying Ea. l|23|l . Ti, Di being determined from numerical fit. The noise rji 
is thus an Ornstein-Ulhenbeck process with properties. 

(?7,(i)) = and {n,{t)7j,{t')) = 6,, — exp ( = 1,2,3....) (24) 

Ti \ Ti J 

Clearly Ti and Di are the correlation time and strength of the color noise variable rji. The c-number noise f{t) due 
to heat bath is therefore given by 

n 

fit) = E (25) 

i=l 

Having obtained the scheme for generation of c-number noise f(t) we now proceed to solve the stochastic differential 
equations. 

In order to solve the c-number Langevin equation Ea. (|12|l we may write it in the equivalent forr" ^^'^^''^^ 



q ^ p 

p = -v'iq) + Qit) + Y,mit) + z 

i 

z = -T^~- (26) 

Tc Tc 
Ti Ti 

The integration of the above set of equationsSSiSS is carried out using the second order Heun's algorithm. A very 
small time step size, 0.001, has been used. 

The above set of equations differ from the corresponding classical equations in two ways. First, the noise correlation 
of c-number heat bath variables f{t) are quantum mechanical in character which is reflected in Di and Ti. Second, 
the knowledge of Q{t) requires the quantum correction equations given in the Appendix A which provides quantum 
dispersion about the quantum mechanical mean values q and p of the system. It is thus essential to take care of 
these contributions. A simple way to achieve this is to solve the equations Ea. (|Al|l to Ea. l|A3p in the Appendix A 
for {5q^(t)) by starting with N-particles (say, around 5000) all of them above the barrier top at g = 0, half with a 

positive velocity distributed according to velocity distribution's p exp (^ 2hcuo<Jk:+i/'^) ) ^^"^ other half with the same 
distribution with negative velocities, initial values of dispersion being set as {Sq'^{t))t=o = 0.5, {6q Sp+6p Sq)t=o = 1.0, 
{Sp^(t))t=o — 0.5 and with others set as zero. The width of the distribution is the same as that for Eq.® where utj is 
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replaced by uq corresponding to the reactant harmonic well which is in equilibrium with the bath. We take the time 
averaged contribution of the quantum corrections upto l/F time for each trajectory to solve the Langevin equation 
Ea. H26|l for N-particles. The time dependent transmission coefficient is calculated from these sets of simulated data 
by calculating 



<t) ^^^ILR (27) 

where N+{t) and N^{t) are the particles that started with positive velocities and negative velocities, respectively 
and at time t are in or over the right hand well (i.e. the particles for which the quantum mean value q{t) > 0). 



IV. RESULTS AND DISCUSSIONS 



We now consider the potential of the form V{q) ^ a — h , where a and h are the two parameters. The other 
three input parameters for our calculations are temperature T, strength of noise correlation F and correlation time 
Tc- For the present purpose we fix a = 0.001 and b ~ 0.5 for the entire calculation except those for Fig. 3 and Fig. 8. 
In order to ensure the stability of the algorithm we have kept Ai/rc << 1 where Ai is the integration time step size. 
Fig. 2 exhibits the temporal variation of classical transmission coefficient (dotted line) and transmission coefficient 
calculated by present method (solid line) for two typical different parameter regimes (a) F = 3.0 and (b) F = 5.0 for 
kT = 0.5, Tc = 3.0 to illustrate the differential behaviour of the classical2& and quantum effects. In the both cases one 
observes significant increase in transmission coefficient due to quantum contribution over and above the magnitude of 
classical transmission coefficient. In order to extract out the contribution of quantum correction due to nonlinearity 
of the system potential Q{t) we exhibit in Fig. 3 the time dependent transmission coefficient for the parameter set 
a = 0.005, h = 0.5, kT = 0.025, F = 1.0 and Tc = 1.0 with (dotted line) and without (solid fine) second order quantum 
corrections. It is observed that the quantum corrections due to nonlinearity affects the stationary values more than 
the transient ones. 

In order to check the workability of the method we now compare in Fig. 4 the temporal variation of numerical 
(dotted line) and analyticaliS (solid line) transmission coefficient at kT = 1.0 for two different parameter regimes 
characteristic of adiabatic regime (a) F — 2.0, Tc — 5.0 (b) F = 4.0, Tc — 8.0 and caging regime (c) F — 90.0, Tc — 10.0, 
the two regimes being differentiated according to lo^, F and Tc; adiabatic (w^ — T/tc > 0); caging ((w^ — F/tc < 0), 
ujb refers to the barrier frequency). Wc observe that while in the adiabatic regime the agreement is excellent the 
numerical transmission coefficients tend to settle down around zero relatively earlier compared to analytical one in 
the caging regime; the numerically calculated phase of oscillation, however, corresponds correctly to its analytical 
counterpart. Keeping in view of the fact that the analytical results are based on a phase space function approach 
which is independent of the method of numerical simulation, we believe that the agreement is quite satisfactory in 
both adiabatic and caging regimes which ensures the validity of the numerical procedure as followed in the present 
treatment. 

In Fig. 5 we show the typical simulation results over a range of dissipation parameters F from 0.01 to 8.0 for Tc = 1.0 
and kT = 1.0. For a low value of F, e.g., 0.01 simulation result remains at a high value for some period to drop rather 
suddenly and to settle finally in an oscillatory way to a low asymptotic value. The behaviour is almost qualitatively 
same as its classical counterpart as shown by Sancho, Romero and Lindenberg2&. The long temporal oscillation of 
transmission coefficient calculated by the c-number method is typical for very low dissipation regime. As dissipation 
increases to 0.5 the temporal variation of K{t) becomes monotonic and K{t) settles down much earlier at a much higher 
asymptotic value. For F = 1.0 to 8.0 the transmission coefficient decreases rather significantly. The oscillations at 
short times for F = 5.0 and 8.0 are due to usual transient recrossings. 

The variation of asymptotic transmission coefficient with dissipation parameter F in Fig. 5 illustrates the well known 
turnover phenomenon. In Figs. 6-8 we analyze this aspect in greater detail. To this end we show in Fig. 6 how the 
asymptotic transmission coefficient as a function of dissipation constant F calculated (classical) earlier by Lindenberg 
et aP^ follows closely to that of ours at relatively high temperature kT = 3.0 for Tc = 3.0. The agreement between the 
two allows us to have a numerical check on the method and to ensure the validity of the result in the classical limit. 
The turnover problem was investigated earlier by Melnikov et aP^ and PoUak et alS^ to obtain correct analytical 
solution bridging the spatial diffusion and the energy diffusion regimes and that goes over to the limiting behaviour 
at high and low dissipation. To explore the quantum effects we study the turnover behaviour at various temperatures 
down to absolute zero as shown in Fig. 7 for Tc = 3.0. As the temperature is lowered the maximum at which the 
turnover occurs shifts to the left and the damping regime that corresponds to classical energy diffusion, i.e., in the 
low damping regime becomes exponentially small as one approaches to absolute zero. It is fully consistent with the 
earlier observation^ made on this issue. 
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To check the vahdity of the numerical results on quantum turnover in a more quantitative way we have further 
compared in Fig. 8 our results (solid lines) with full quantum results (dotted lines) of Topaler and Makriii (Fig. 
9a-b of Ref.l2) for two different scaled temperatures kT = 1.744(200iv:) and kT = 2.617(300X) for a double well 
potential with a — 0.0024 and h = 0.5 in the Ohmic regime. It is immediately apparent from Fig. 8 that the 
transmission coefficients calculated by simulation of the c-number Langevin equation with fourth order quantum 
corrections compared satisfactorily with the full quantum results based on path integral Monte Carlo method of 
Topaler and Makri^^. 

In Fig. 9 we plot the temporal behaviour of transmission coefficient K[t) for different temperatures for F = 1.0 
and Tc — 1.0. It is observed that K{t) quickly settles after a fast fall and the asymptotic values of the transmission 
coefficient converge to a temperature independent value as the temperature is increased from kT = 0.0 to kT = 5.0. 
Fig. 10 exhibits this variation of asymptotic value of transmission coefficient (dotted line) explicitly as a function of 
temperature and a comparison with analyticaliS results (solid line) for F — 2.0 and — 5.0. The agreement is found 
to be excellent. In order to analyze the temperature dependence of the transmission coefficient at low temperature 
and compare with earlier results^k^ the numerical result (solid line) is fitted against a function of the form (dotted 
line) Aexp{B /T^), in the inset of Fig. 10, where A and B are fitting constants. The well known T^ enhancement of the 
quantum rate is observed upto a low temperature beyond which the fitting curve tends to diverge rapidly as kT 0. 
This divergence of the rate at very low temperature had been noted earlier in the analytical result of Grabert et ai^. 
In the present calculation, however the transmission coefficient reaches its maximum finite value within unity and its 
validity is retained even in the vacuum limit. Thus in contrast to classical transmission coefficient the temperature 
dependence remains a hallmark of quantum signature of the transmission coefficient. In Fig. 11 we show the temporal 
behaviour of the transmission coefficient for several values of the noise correlation time Tc for kT = 1.0 F = 1.0. The 
asymptotic transmission coefficient increases as expected from theoretical point of view. 

V. CONCLUSION 

The primary aim of this paper is to extend our recent treatment of c-number Langevin equation to calculate 
numerically the time dependent transmission coefficient within the framework of reactive flux formalism. Since the 
quantum dynamics is amenable to a theoretical description in terms of ordinary coupled equations which are classical 
looking in form it is possible to employ the numerical simulation scheme of Straub and Berne in a quite straight forward 
way for calculation of transmission coefficient. There are two special advantages in the procedure. First, since we 
are concerned with the dynamics of quantum mechanical mean values coupled to quantum correction equations, we 
need not calculate higher moments {<f) or (p^) etc in any stage of calculations. This makes the calculation simpler. 
Secondly, the treatment can be readily used to calculate the dynamics even in the vacuum limit, i.e., kT — > 0, where 
it is expected that because of the oscillating nature of the real time propagator in path integral method Monte Carlo 
schemes pose very serious problems from applicational point of view. We have extended the classical simulation 
procedure to a quantum domain taking into consideration of the quantum effects in two different ways. The quantum 
effects enter through the correlation function of the c-number noise variables of the heat bath and furthermore, through 
nonlinearity of the system potential which is entangled with quantum dispersion around the quantum mechanical mean 
values of the system operators. We summarize the main results as follows: 

(i) The present method is a direct extension of classical simulation method of Straub and Berne to quantum domain 
for calculation of transmission coefficient within a c-number version of reactive flux formalism. Although the quantum 
effects due to heat bath can be taken into account in terms of noise correlation expressed in quantum fluctuation- 
dissipation relation, the quantum dispersion around the quantum mean values of the system operators are to be 
calculated order by order. Notwithstanding the latter consideration the method is efficient when compared to full 
quantum mechanical calculation as demonstrated in the present simulation. 

(ii) We have calculated the time dependent transmission coefficient for a double well potential over a wide range of 
friction and temperature and shown that our numerical simulation results on turnover phenomena, low temperature 
enhancement of quantum rate and other features agree satisfactorily with those calculated analytically/otherwise using 
phase space function and other approaches. The differential behaviour of the classical and transmission coefficients 
calculated by the classical and present c-number method has been analyzed in detail. The procedure is equipped to 
deal with arbitrary noise correlation, strength of dissipation and temperature down to vacuum limit, i.e., kT 0. 

(iii) The present approach is independent of path integral approaches and is based on canonical quantization and 
Wigner quasiclassical phase space function and takes into account of the quantum effects upto a significant degree of 
accuracy. This procedure simply depends on the solutions of coupled ordinary differential equations rather than the 
multi-dimensional path integral Monte Carlo techniques and is therefore complementary to these approaches, much 
simple to handle and corresponds more closely to classical procedure. 
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APPENDIX A 

Evolution Equations For Higher-Order Quantum Corrections For Anharmonic Potential 

The equations upto fourth order for quantum corrections (corresponding to the contribution of anharmonicity of 
the potential) with dissipative effects taking into consideration of Lorentzian density of reservoir modes of the form 
Ea. (|20|l . in the limit very small, are listed below. 
Equations for the second order are: 

J^{5f) = {SqSp + SpSq), 

^(Sp^) = -2T(6f) - V"(6qSp + SpSq) - V'"(SqSp6q), (Al) 
at 

^{Sq6p + SpSq) = -T{SqSp + Sp5q)2{6f)-2V"{df)~V'"{Sf), 
Those for the third order are: 

j^iSf) - 3{SqSp6q), 

j^{6f) = -3r{Sf) - 3V"{SpSq6p) + V" (^^{Sf){Sp') - ^{SpSfSp) + 

-{SqSpSq) = -T{SqSp5q) + 2{SpSqSp) - V" (Sf) - — [{Sq') - (dqY) , (A2) 
^{6pSqSp) ^ -2r (SpSqSp) + {6f) - 2V"{6q6p5q) 



dt 



ylll 

— {{6f)(SqSp + Sp6q) - {Sq^Sp + 5pSf)) 



And lastly, the fourth order equations are: 



j^iSq') = 2{df6p + 6pdf), 

= -^T{6p^) - 2V"{5q5f + SfSq) + 2V"'{5q^){5p^), 
^{6f6p + SpSf) = -T{6fSp + 6p6q^)-2V"{6q'^)~3h^ + 6{SpSfSp) 



+ V"'{Sf){Sf), (A3) 
{SqSf + 5f5q) = -3r{dq5f + Sp^Sq) + 2{Sp^) + W" {fi'^ - 2{5p6fSp)) 



dt 

+ W"{5q^){5p5q5p), 
^;{SpdfSp) = -2T{6pdfSp) - V"{dfdp + 6pSf) + {SfSq + SqSp^) 



dt 

+ V'"{6f){SqSpSq)- 

The derivatives of V{q), i.e., V" or V'" etc. in the above expressions are functions of q the dynamics of which is given 
by Eq. fmi . 
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FIGURE CAPTIONS 

Fig.l: Plot of correlation function c{t) vs t as given by Eq.(19) (solid line) and Eq.(21) (dotted line) for the set of 
parameter values mentioned in the text. 

Fig. 2: A comparison of transmission coefficients ( classical(dotted line); c-number method(solid line)) as function 
of time t is plotted for (a) T = 3.0 (b) F = 5.0 for the parameter set mentioned in the text. 

Fig. 3: The numerical transmission coefficient K{t) is plotted against time with (dotted line) and without (solid line) 
quantum correction Qit) for a = 0.005, b = 0.5, F = 1.0, Tc = 1.0 at kT = 0.025 

Fig. 4: Numerical transmission coefficient K{t) is plotted against time t and compared with c-number analyticaliS 
K{t) for adiabatic regime [(a) F = 2.0, Tc = 5.0, (b) F = 4.0, Tc = 8.0] and caging regime [ F = 90.0, Tc = 10.0] for the 
parameter set mentioned in the text. 

Fig. 5: Transmission coefficient, K{t) is plotted against time for different values of F for the parameter set mentioned 
in the text. 

Fig. 6: A comparison of the turnover (plot of asymptotic k vs F) calculated by Sancho, Romero and Lindenberg2& 
(circle) with that by the present method (square) for the parameter set mentioned in the text. 

Fig. 7: Kramers' turnover (plot of asymptotic k vs F) for different temperatures for the parameters set mentioned 
in the text for kT — 0.0 (downtriangle), kT = 0.5 (circle), kT = 1.0 (square) and kT = 3.0 (uptriangle). 

Fig. 8: A comparison of the Kramers' turnover (plot of assymptotic k vs F) for two temperatures (a) kT — 
2.617(300/4:) and (b) kT = 1.744(200/4:) between the present result (sohd line) and full quantum result (dotted 
line) of Topaler and Makri (Fig.9a-b of ref.l2) for the double well potential with a = 0.0024 and b = 0.5 in the Ohmic 
regime. 

Fig. 9: Transmission coefficient K{t) is plotted against time t for different values of temperature, kT — 0.0 (solid 
line), kT = 0.5 (dash dot dot line), kT = 1.0 (dash dot line), kT = 2.0 (dashed line) and kT = 5.0 (dotted line) for 
the set of parameters mentioned in the text. 

Fig. 10: Asymptotic transmission coefficient is plotted against temperature (dotted line: analytical^^, solid line: 
numerical) for the parameter set mentioned in the text. Inset: the same numerical curve (solid line) is plotted against 
a fitted curve (dotted) to exhibit T^ enhancement of rate at low temperature. 

Fig. 11: Transmission coefficient K{t) is plotted against time for different values of Tc for the parameter set mentioned 
in the text. 
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